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Abstract:  A  depth-averaged  two-dimensional  nonuniform  sediment  transport  model 
is  applied  to  the  beaches  adjacent  to  Grays  Harbor,  WA,  USA  to  test  the  model  skill  in 
predicting  nearshore  morphology  change.  The  model  considers  bed  material  hiding, 
exposure,  sorting,  stratification,  bed-slope  effects,  avalanching,  non-erodible  bed 
surfaces,  and  transport  due  to  asymmetrical  waves,  Stokes  drift,  roller  and  undertow. 

The  sediment  transport,  bed  change  and  sorting  equations  are  solved  simultaneously 
and  implicitly  at  the  same  time  step  as  the  hydrodynamics.  The  model  is  able  to  capture 
the  onshore  migration  of  the  offshore  bar  and  filling  of  the  trough  but  has  difficulty  in 
the  foreshore  region  where  swash  zone  processes  are  neglected.  The  calculated 
nearshore  water  depths  agree  with  measurements  with  an  average  Brier  Skill  Scores  of 
0.3  and  bed  changes  with  an  average  correlation  coefficient  of  0.53. 


Introduction 

The  grain  size  distribution  of  coastal  sediments  is  a  direct  consequence  of  the 
sediment  sources  and  hydrodynamic  conditions  that  exist  in  each  specific 
environment.  Mason  and  Folk  (1958)  showed  that  sediments  in  different  coastal 
environments  can  be  distinguished  by  their  statistical  parameters  such  as  mean, 
standard  deviation  (sorting)  and  skewness.  For  most  beaches,  coarser  sediment  is 
generally  found  in  the  swash  zone  and  the  wave  breaker  line  while  finer  sediment  is 
found  in  the  trough,  landward  of  the  breaker  line  (e.g.  Mason  and  Folk,  1958;  Wang 
et  al.,  1998).  As  reported  in  Mason  and  Folk  (1958),  in  most  beaches  the  sediments 
are  negatively  skewed  (coarser  grained)  because  the  fine  sediments  are  winnowed 
away  by  breaking  waves.  In  general,  finer  sediments  are  found  in  areas  of  flow 
deceleration  and  vice-versa  (e.g.  Black  et  al.,  1989).  Many  coastal  inlets  have  shell- 
lag  deposits  in  the  inlet  throat  which  armors  the  and  prevents  excessive  erosion.  Bed 
material  textural  changes  can  also  be  related  to  storm  events,  seasonal  climatic 
changes  and  long-term  depositional  and  erosional  trends  due  to  changes  in  the 
amount  or  properties  of  the  sediment  source(s).  For  example,  Wang  et  al.  (1998) 


*  All  figures  are  in  color  on  the  Coastal  Sediments  Proceedings  DVD. 
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found  that  storms  and  the  resulting  offshore  migration  of  the  bar  could  leave  a  layer 
of  coarser  lag  deposit  where  fine  deposits  would  normally  be  found. 

The  Coastal  Modeling  System  (CMS)  developed  by  the  U.S.  Army  Corps  of 
Engineers  (Buttolph  et  al.,  2006;  Lin  et  ah,  2008)  is  an  integrated  wave, 
hydrodynamic,  sediment  transport  and  morphology  change  modeling  system.  The 
CMS  has  previously  been  validated  at  Grays  Harbor  for  hydrodynamics  and  waves 
in  Lin  et  al.  (2008)  and  Wu  et  al.  (2010, 2011).  In  this  paper,  a  depth-averaged  two- 
dimensional  (2DH)  nonuniform  total  load  sediment  transport  model  is  added  to  the 
CMS  and  applied  to  Grays  Harbor,  WA,  USA  in  order  to  test  the  model 
performance  in  predicting  nearshore  morphology  change  and  distribution  of 
nearshore  sediments.  The  model  uses  a  multi-fraction  method  in  which  the  sediment 
mixture  is  divided  into  discrete  fractions  and  the  total  sediment  transport  rate  is 
obtained  as  the  sum  of  the  individual  fractional  transport  rates.  The  bed  is  also 
divided  into  discrete  layers  and  the  fractional  composition  of  each  sediment  layer  is 
simulated  in  time.  The  model  also  includes  transport  over  nonerodible  bed  surfaces, 
avalanching,  and  cross-shore  transport  due  to  wave  asymmetry,  roller,  and 
undertow. 


Nonuniform  Sediment  Transport  Model 

The  single-sized  sediment  transport  model  described  in  Sanchez  and  Wu  (2011)  is 
extended  to  multiple- sized  nonuniform  sediments.  In  this  model,  the  sediment 
transport  is  separated  into  current-  and  wave-related  transports.  The  transport  due  to 
currents  includes  the  stirring  effect  of  waves;  and  the  wave-related  transport 
includes  the  transport  due  to  asymmetric  oscillatory  wave  motion  and  also  steady 
contributions  by  Stokes  drift,  surface  roller,  and  undertow.  The  current-related  bed 
and  suspended  transports  are  combined  into  a  single  total-load  transport  equation, 
thus  reducing  the  computational  costs  and  simplifying  the  bed  change  computation. 
The  2DH  transport  equation  for  the  current-related  total  load  is 


dt 
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for  7  =  1, 2;  k  =  l,2,...,N ,  where  N  is  the  number  of  sediment  size  classes,  t  is  time, 
h  is  the  total  water  depth,  Xj  is  the  Cartesian  coordinate  in  the  jth  direction,  Uj  is  the 
depth-averaged  current  velocity,  Ctk  is  the  depth-averaged  total-load  sediment 
concentration,  f^tk  is  the  total-load  correction  factor,  is  the  fraction  of  suspended 
load  in  total  load  (equal  to  1  for  fine-grained  sediments),  is  the  sediment  mixing 
coefficient,  is  the  total-load  adaptation  coefficient  and  0)^^.  is  the  sediment  fall 
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velocity.  The  equilibrium  concentration  is  calculated  as  where  is 

fraction  of  the  k'^  sediment  size  class  in  the  top-most  bed  layer  and  is  the 
potential  equilibrium  sediment  concentration  calculated  from  empirical  formulas. 
The  hiding  and  exposure  of  each  bed  material  size  class  are  considered  by 
modifying  the  critical  shields  parameter  using  the  method  of  Wu  et  al.  (2000). 

The  correction  factor  f^tk  is  the  ratio  of  the  depth-averaged  sediment  and  flow 
velocities  and  accounts  for  the  vertical  nonuniform  profiles  of  sediment 
concentration  and  current  velocity.  In  a  combined  bed  load  and  suspended  load 
model,  the  correction  factor  is  given  by  =  ^l\jsk  I  f^sk  +  (1  “  ^sk  )U  /  ]  where 

Ubk  is  the  bed  load  velocity  and  fitk  is  the  suspended  load  correction  factor  and  is 
defined  as  the  ratio  of  the  depth-averaged  sediment  and  flow  velocities.  Because 
most  of  the  sediment  is  transported  near  the  bed,  both  the  total  and  suspended  load 
correction  factors  are  usually  less  than  1  and  typically  in  the  range  of  0.3  and  0.7, 
respectively.  By  assuming  logarithmic  current  velocity  and  exponential  suspended 
sediment  concentration  profiles,  an  explicit  expression  for  the  suspended  load 
correction  factor  may  be  obtained  as 


E,  (IPA)  -  E,  (^)  +  ln(.4  /  Z)  t*''  -  ln(l  /  Z)  e'** 
e"**"*  [ln(l  /  Z)  - 1]  [l  -  ] 

where  (f>-  (O^h!  e  ,  A  -  a!  h  ,  Z  =  Zq!  h  ,  £  is  the  vertical  mixing  coefficient,  and 
Ej  is  the  exponential  integral.  The  equation  can  be  further  simplified  by  assuming 
that  the  reference  height  is  proportional  to  the  roughness  height  (e.g.  a  =  SOzq  ),  so 
that  (Z,^) . 


The  fractional  bed  change  is  calculated  as 


(i-p™) 


^l^sk  {^Ik*  ^tk 


dQ„k.j  ^  d 

dxj  dxj 


K] 

dxj^ 
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where  ^  is  the  still  water  depth,  p'^  is  the  bed  porosity,  Q^k  is  the  sediment 
transport  due  to  wave  asymmetry,  Stokes  Drift,  roller  and  undertow,  D,  is  a  bed- 
slope  coefficient,  and  Qbk  =  hUCtk{l-rsk)  is  the  bed  load.  The  total  bed  change  is 
calculated  as  the  sum  of  Eq.  (2)  for  all  size  classes.  The  advantage  of  including  Q^k 
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in  the  bed  change  equation  instead  of  the  transport  equation  is  that  it  simplifies  the 
calculation  and  has  been  found  to  improve  model  stability.  In  addition,  it  is  not 
straightforward  to  include  the  wave-related  transport  in  Eq.  (1),  since  modifications 
would  need  to  be  made  to  the  total  load  correction  factor  and  adaptation  coefficient. 
For  simplicity,  it  is  assumed  that  the  onshore  and  offshore  components  are  on  the 
same  axis  as  the  waves.  Details  on  are  left  for  a  later  publication. 

The  bed  material  above  the  erodible  layer  is  divided  into  multiple  layers,  and  the 
sorting  of  sediments  is  calculated  using  the  mixing  layer  concept.  The  mixing  layer 
is  the  top-most  (first)  layer  of  the  bed.  The  temporal  variation  of  the  bed-material 
composition  in  each  layer  is  calculated  using  a  method  similar  to  Wu  (2004).  The 
sediment  transport,  bed  change,  and  bed  gradation  are  solved  simultaneously 
(coupled),  but  are  decoupled  from  the  flow  calculation  at  the  time  step  level.  To 
illustrate  the  bed  layering  process.  Figure  1  shows  an  example  of  the  temporal 
evolution  of  7  bed  layers  during  erosional  and  depositional  regimes.  Details  on  the 
mixing  layer  thickness  calculation  and  the  bed  layering  algorithm  are  left  for  a 
subsequent  publication. 


Fig.  1 .  Schematic  showing  an  example  bed  layer  evolution.  Colors  indicate  layer  number. 


Field  Study 

Grays  Harbor  inlet,  WA  is  located  on  the  southwest  Washington  coast,  USA,  at  the 
mouth  of  the  Chehalis  River.  Between  May  and  July  of  2001,  the  U.S.  Geological 
Survey  instrumented  6  tripods  and  collected  time  series  of  wave  height,  water 
surface  elevation,  near-bottom  current  velocity,  and  sediment  concentration  proxies 
(Fanderman  et  ah,  2004).  Weekly  topographic  maps  and  monthly  bathymetric 
surveys  along  transects  spaced  50-200  m  apart  were  collected  (see  Figure  2).  In 
addition  grab  sample  of  surface  sediment  were  collected  at  several  locations.  Figure 
1  shows  the  location  of  the  observation  stations  and  monthly  nearshore  bathymetric 
profiles. 
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Fig.  2.  Map  of  Grays  Harbor  inlet,  WA  showing  the  location  of  the  nearshore  bathymetric  transects. 

Model  Setup 

The  first  half  of  the  field  deployment  between  May  6-30  of  2001  was  simulated. 
The  simulation  period  was  characterized  by  relatively  calm  conditions,  with  a  few 
spring  storms  with  significant  wave  heights  on  the  order  of  3  m  (see  Figure  3). The 
spectral  wave  transformation  model  CMS-Wave  was  run  on  a  -200,000  cell 
Cartesian  grid  with  varying  grid  resolution  from  15-120  m  (see  Figure  4).  The 
waves  were  forced  with  spectral  wave  information  from  the  Coastal  Data 
Information  Program  (CDIP)  buoy  No.  03601  located  southwest  of  the  inlet  at  a 
depth  of  42  m  (see  Figure  2). 


Fig.  3.  Time  series  of  significant  wave  height  (top),  peak  wave  period  (middle),  and  incident  wave  angle 
with  respect  to  shore  line  (bottom)  during  May  2001. 


The  CMS-Flow  was  forced  with  a  water  level  time  series  from  Westport  Harbor 
with  a  negative  30  min  phase  lag  correction  which  was  obtained  by  comparing  with 
the  stations  deployed  during  the  field  study  (see  Figure  5).  Winds  were  interpolated 
from  the  Blended  Sea  Winds  product  of  the  National  Climatic  Data  Center  (Zhang 
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et  al.  2006).  The  Manning’s  coefficient  was  calibrated  in  previous  studies  as  0.018 
over  the  whole  domain  except  on  the  rock  structures  where  a  value  of  0. 1  was  used. 
A  flux  boundary  condition  was  applied  at  the  Chehalis  River  which  was  obtained 
from  the  USGS.  The  CMS -Flow  ~5 5, 000-cell  quadtree  grid  is  shown  in  Figure  6 
and  has  six  levels  of  refinement  from  20-640  m.  A  variable  time  step  was  used  with 
a  maximum  value  of  10  min.  The  sediment  transport  and  bed  change  were 
calculated  at  every  hydrodynamic  time  step. 


Fig.  4.  Nonuniform  Cartesian  grid  used  for  CMS-Wave. 


A  ramp  or  spin  up  period  of  5  days  was  implemented  based  on  previous  hydrodynamic 
studies  at  Grays  Harbor  so  the  start  of  the  simulation  was  May  1,  2001.  Waves  were 
calculated  at  a  constant  2  hr  interval  (steering  interval).  The  significant  wave  height, 
peak  wave  period,  wave  unit  vectors,  and  wave  dissipation  were  linearly  interpolated 
to  the  flow  grid  every  steering  interval  and  then  linearly  interpolated  in  time  at  every 
hydrodynamic  time  step.  Wave  variables  such  as  wave  length  and  bottom  orbital 
velocities  were  updated  every  time  step  for  wave-current  interaction.  When  using  such 
a  large  steering  interval,  it  is  important  to  consider  how  the  water  levels,  current 
velocities  and  bed  elevations,  which  are  passed  from  the  flow  to  the  wave  model,  are 
estimated.  For  this  application,  and  for  most  open  coast  applications,  the  nearshore 
waves  are  most  sensitive  to  variations  in  water  levels  and  not  currents.  Therefore, 
improved  results  can  be  obtained  by  predicting  the  water  levels  at  the  wave  model 
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time  step  based  on  a  decomposition  of  the  water  levels  into  spatially  constant  and 
variable  components.  The  spatially  constant  component  is  assumed  to  be  equal  to  the 
tidal  water  surface  elevation  and  the  spatially  variable  component  which  includes  wind 
and  wave  setup  is  estimated  based  on  the  last  flow  time  step.  The  currents  and  bed 
elevations  which  are  passed  from  the  wave  to  flow  grid  are  simply  set  to  the  last  time 
step  value.  Other  types  of  prediction  methods  could  be  used;  however,  the  approach 
described  above  has  been  found  to  be  sufficient  for  most  applications  and  is  simple  to 
calculate.  After  each  wave  run,  a  surface  roller  model  is  also  calculated  on  the  wave 
grid  and  the  roller  stresses  are  added  to  the  wave  stresses  before  interpolating  on  to  the 
flow  grid.  Even  though  CMS -Flow  and  CMS -Wave  use  different  grids,  the  two 
models  are  in  a  single  code  which  facilitates  the  model  coupling  and  speeds  up  the 
computation  by  avoiding  communication  files,  variable  allocation  and  model 
initialization  at  every  steering  interval. 


Fig.  5.  Water  levels  (top)  and  wind  velocities  and  direction  (bottom)  during  May  of  2001. 


Fig.  6.  Quadtree  grid  used  for  CMS-Flow. 
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The  initial  bed  material  composition  was  specified  by  a  spatially  variable  median 
grain  size  and  constant  geometric  standard  deviation  of  1.3  mm  based  on 

field  measurements.  The  initial  fractional  composition  at  each  cell  was  assumed  to 
be  constant  in  depth  and  have  a  log-normal  distribution,  and  represented  by  six  size 
classes  with  characteristic  diameters  of  0. 1 , 0. 126, 0.16, 0.2, 0.25,  and  0.3 1  mm.  An 
example  of  the  initial  grain  size  distribution  is  shown  in  Figure  7.  Ten  bed  layers 
were  used  with  an  initial  thickness  of  0.5  m  each.  The  Lund-CIRP  transport 
formulas  were  used  to  estimate  the  transport  capacity  (Camenen  and  Larson,  2007). 
The  total-load  adaptation  coefficient  is  calculated  as  =Uhl  (L^O)^ )  where  is 
the  total-load  adaptation  length.  Here  =  (1-  -I-  ,  where  and  are 

the  bed-  and  suspended-load  adaptation  lengths,  respectively.  The  bed-load 
adaptation  length  is  set  to  10  m,  and  the  suspended-load  adaptation  length  is 
calculated  as  =Uhl  where  the  suspended-load  adaptation  coefficient 

is  set  to  0.5.  A  constant  bed  porosity  of  0.3  was  used  in  the  simulation.  The  bed  load 
velocity  is  calculated  using  Van  Rijn’s  (1984)  formula  with  recalibrated  coefficients 
by  Wu  et  al.  (2006). 


Grain  size,  mm 

Fig.  7.  Example  log-normal  grain  size  distribution  (  d^Q  =  0.16  mm,  =  1.3  mm). 


Results  and  Discussion 

Calculations  were  performed  on  a  desktop  PC  and  the  31  days  simulation  was 
completed  in  approximately  10  hrs.  A  comparison  of  the  measured  and  computed 
bed  changes  between  May  6  and  30  of  2001  is  shown  in  Figure  8.  Selected  regions 
of  interest  are  encompassed  by  black  lines  in  order  to  help  visually  compare  the  bed 
changes.  In  general,  the  results  show  many  common  features  and  similar  erosion 
and  deposition  patterns.  More  specifically,  the  bed  change  is  characterized  by  the 
erosion  of  the  outer  bar,  deposition  in  the  inner  bar  face  and  outer  trough,  and 
erosion  of  the  inner  trough  face.  There  is  a  region  extending  approximately  1  km 
from  the  northern  jetty,  where  the  bed  changes  are  noticeably  different  from  those 
further  to  the  north.  This  region  is  interpreted  as  being  strongly  influenced  by  the 
presence  of  the  inlet,  ebb  shoal  and  northern  jetty.  As  an  example.  Figure  9  shows  a 
snap  shot  of  the  current  velocities  on  May  14,  2001  during  an  ebb  tide  and 
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southeasterly  swell  in  which  there  is  a  reverseal  in  the  longshore  current  in  the 
transition  region  due  to  the  combined  effect  of  refraction  of  waves  over  the  ebb 
shoal,  deflection  of  the  ebb  current  by  the  northern  jetty,  and  formation  of  the 
northern  ebb  jet  gyre.  This  emphasizes  the  importance  of  accurately  capturing  the 
inlet  dynamics  in  modeling  adjacent  beaches.  Interestingly,  both  the  measurements 
and  model  results  show  small  (200-300  m  in  length)  inner  bars  form  adjacent  to  the 
trough,  which  appear  to  occur  at  regular  400-500  m  intervals. 

The  computed  bed  changes  in  the  foreshore  region  (beach  face)  are  relatively  small 
compared  to  the  measurements  due  to  the  lack  of  swash  zone  processes  in  the 
present  version  of  CMS.  Swash  zone  processes  enhance  transport  in  the  surf  zone 
by  increasing  the  current  velocities,  transport  rates  and  mixing  at  the  shoreline. 
Large  component  of  longshore  sediment  transport  occurs  in  the  swash  zone  and 
without  these  processes,  morphodynamic  models  will  tend  to  underestimate 
longshore  transport  rates  and  bed  change  in  the  foreshore  region..  Walstra  et  al. 
(2005)  simulated  the  bed  change  at  transects  9  and  20  using  a  two-dimensional 
vertical  (2DV)  profile  evolution  model  and  were  able  to  predict  the  onshore 
migration  of  the  bar,  but  also  found  that  the  model  performance  deteriorates  in  the 
foreshore  region. 

The  measured  and  computed  water  depths  and  bed  changes  for  Transects  1  and  9 
are  shown  in  Figure  10.  As  observed  in  Figure  8,  most  of  bed  changes  occurred 
from  the  nearshore  bar  to  the  outer  beach  face.  The  model  was  able  to  accurately 
predict  an  onshore  bar  migration  although  it  underestimated  the  nearshore  bar 
height  which  is  also  observed  in  Figure  8.  In  order  to  evaluate  the  model 
performance  in  predicting  the  nearshore  bathymetry,  the  Brier  Skill  Score  BSS  was 
applied  to  the  water  depths  and  the  correlation  coefficient  to  the  bed  change. 
Other  goodness  of  fit  parameters  were  also  calculated  and  showed  similar  patterns, 
for  simplicity  only  the  previously  mentioned  parameters  are  shown  in  Figure  10. 
The  goodness  of  fit  statistics  show  a  wide  range  of  values. 

The  measured  bed  change  shows  a  larger  variation  as  compared  to  the  model 
results,  indicating  that  morphology  change  is  sensitive  to  longshore  variations  in 
forcing,  initial  bathymetry  or  3D  processes  such  as  rip  currents.  As  discussed  by 
Walstra  et  al.  (2005),  the  model  results  indicate  that  the  waves  and  currents  do  in 
fact  vary  over  the  spatial  scales  (10-100  m)  of  the  observed  morphological 
variations.  As  an  example.  Figure  12  shows  a  snap  shot  of  the  current  velocities 
during  an  ebb  tide  and  easterly  wave  event  (approximately  normal  to  shore)  on  May 
6,  2001  which  shows  complex  nearshore  currents,  several  reversals  in  the  direction 
of  the  longshore  current  and  the  influence  of  the  ebb  jet  gyre. 
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Fig.  9.  Snapshot  of  current  magnitudes  during  an  ebb  tide  on  May  14,  2001. 


Although  the  field  measurements  show  evidence  of  stratification,  no  measurements 
were  conducted  in  the  surf  zone  where  most  of  the  bed  change  occurred.  It  is 
expected  that  due  to  the  strong  mixing  in  the  surf  zone,  the  2DH  model  is  a 
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reasonable  approximation.  In  addition,  previous  studies  at  Grays  Harbor  (Wu  et  al. 
2010,  2011)  have  found  good  agreement  with  hydrodynamic  measurements  and 
sediment  transport  patterns  using  a  2DH  approach. 


Fig.  10.  Measured  and  computed  water  depths  (top)  and  bed  changes  (bottom)  for  Transect  9. 


Fig.  11.  Brier  Skill  Score  for  water  depths  and  correlation  coefficient  for  computed  bed  changes  at 
selected  Transects. 


Fig.  12.  Snap  shot  of  current  magnitudes  during  an  ebb  tide  on  May  6,  2001. 
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The  computed  median  grain  size  on  May  30,  2001  is  shown  in  Figure  13. 
Qualitatively,  the  results  agree  well  with  field  measurements  and  typical  findings  for 
most  inlets  and  beaches.  Coarser  sediments  are  found  in  the  beach  face  and  breaker 
line  (offshore  bar)  and  finer  sediments  are  found  in  the  trough  and  offshore  of  the 
surf  zone.  In  addition,  coarser  sediments  are  found  in  the  inlet  entrance  and  finer 
sediments  are  found  on  the  periphery  of  the  ebb  shoal.  In  addition,  it  is  noted  that 
the  area  around  the  jetties  highly  armored  due  to  the  strong  currents  and  large  waves 
present,  which  is  also  observed  in  field. 


Fig.  13.  Distribution  of  median  grain  size  calculated  after  the  25-day  simulation. 


Conclusion 

A  depth-averaged  nonuniform  sediment  transport  model  has  been  developed  and 
applied  to  Grays  Harbor,  WA.  Nearshore  measurements  of  bathymetry  were  used  to 
validate  the  model  during  the  period  of  May  6-30, 2001.  Goodness  of  fit  statistics  of 
water  depths  and  bed  changes,  indicate  generally  reasonable  to  good  model 
performance  although  the  model  skill  varied  significantly,  especially  on  the  beach 
face  where  swash  zone  processes  are  likely  important  and  are  not  presently 
represented  in  the  model.  The  measured  bed  change  shows  larger  degree  of 
variability  as  compared  to  model  results  indicating  that  nearshore  morphology  is 
sensitive  to  longshore  variations  in  forcing  and  cross-shore  processes  which  are 
difficult  to  resolve.  Results  also  show  that  there  is  a  region  adjacent  to  the  north 
jetty  (transition  zone)  which  is  strongly  influenced  by  the  presence  of  the  inlet  due 
to  wave  refraction  over  the  ebb-tidal  delta,  ebb  and  flood  currents  including 
detached  eddies,  and  the  presence  of  the  north  jetty. 
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